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By analyzing preliminary experimental measurements of charge-balance functions from the STAR 
Collaboration at the Relativistic-Heavy-Ion Collider (RHIC), it is found that pictures where bal¬ 
ancing charges are produced in a single surge, and therefore separated by a single length scale, are 
inconsistent with data. In contrast, a model that assumes two surges, one associated with the for¬ 
mation of a thermalized quark-gluon plasma and a second associated with hadronization, provides 
a far superior reproduction of the data. A statistical analysis of the model comparison finds that 
the two-surge model best reproduces the data if the charge production from the first surge is similar 
to expectations for equilibrated matter taken from lattice gauge theory. The charges created in the 
first surge appear to separate by approximately one unit of spatial rapidity before emission, while 
charges from the second wave appear to have separated by approximately a half unit or less. 
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I. INTRODUCTION AND THEORY 


A principal goal of colliding heavy-ions at high energy is to verify whether one can create small drops of equilibrated 
quark gluon plasma in the laboratory. Kinetic equilibration is not surprising given the high collision rates and 
is validated by thermal features of the data, especially observables related to collective flow [1 and jet quenching 
[2]. Chemical equilibrium is more difficult to justify. The final-state abundance of hadrons suggests that chemical 
equilibrium was lost just after hadronization when the temperature was near 165 MeV [3]. In contrast, there is scant 
experimental evidence that chemical equilibrium was maintained during the period when T > 200 MeV when the 
matter is expected to be in a phase of strongly interacting quark gluon plasma (QGP). That chemistry is remarkable. 
Counting spins, colors and flavors, there are 36 light degrees of freedom from the up, down and strange quarks, and 
an additional 16 from gluons. Thus, approximately 52 strongly interacting particles should inhabit a volume on the 
order of one thermal wavelength cubed, ^ {hc/Tf. 

At high temperature and zero baryon density, the chemical make-up of the QGP cannot be quantified by counting 
quarks or gluons because they tend to be off-shell or virtual, so their number is not a well-defined observable. However, 
even though the average charge within a volume V is zero, the fluctuations of the charge characterize the degrees of 
freedom. If one considers the three-by-three fluctuation tensor, 


Xab 


QaQb 

v 


(i) 


with a and b referring to the up, down and strange charge, one can gain insight into the chemistry. If up, down and 
strange quarks are good quasi-particles with no inter-quark correlations (a gaseous state), 

xi? GP) = (n-a + ria)Sab , (2) 

where n a is the density of quarks of the given flavor, up, down or strange. Correlations between quarks, such as the 
two being in the same hadron, alter the expression. For a hadron gas the correlations are 

xib ad) = ^2n a q aa q ab , (3) 

a 

where n a is the density of hadron species a which has a charge q aa . As an example, protons contribute a factor of four 
times their density to Xuu because the “up” charge of a proton is two. Even for strongly interacting and correlated 
systems Xab is a well-defined observable, because charge is conserved by the strong interaction. From equation [3] one 
can see that hadronic resonances induce off-diagonal elements to Xab 0 • For example, the K + contributes negatively 
to Xus and a A hyperon gives a positive contribution to Xus • 

For a massless gas of quarks and gluons, the number of quarks within a fluid element of fixed entropy stays constant 
during an isentropic expansion because both the number densities and the entropy density scale as T 3 . Therefore, 
in an isentropic expansion of a massless parton gas the number of quarks within the fluid element stays relatively 
constant. Lattice calculations show that this property remains reasonably preserved even for a strongly interacting 
system. The ratio Xab/s is illustrated in Fig. [l] where it changes only at the 10% level once T > 225 MeV. Below 
that temperature, the system is hadronizing and x is strongly temperature dependent. Due to entropy conservation, 
the number of hadrons just below T c is marginally lower than the number of quarks above T c , so given that each 
hadron has two or three quarks, copious quark production, through string breaking or resonance decays, accompanies 
hadronization. As shown in Fig. [l] some of the elements of the fluctuation tensor change significantly, especially 
Xuu/s « Xdd/s which nearly doubles. 

Unfortunately, charge fluctuations, Xab(T), are not directly accessible from experiment. If one could measure the 
charge fluctuation within a small volume, a few fm 3 , one would expect the fluctuation to approximate that of an 
equilibrated system as long as the equilibrated system does not have correlations beyond that scale. However, if one 
considers the net charge in the entire system, it does not fluctuate due to the fact that, unlike the assumptions for a 
grand canonical ensemble, the net charge is fixed and X ^ 0. In contrast, the charge correlation 


gabi^v) = (Pa(O)pb(AT?)) , (4) 

is sensitive to both the equilibrated charge fluctuation and accounts for the conservation of net charge. Here, p a (jj) * s 
the charge density per unit of spatial rapidity and the brackets denote averages over events. Translational invariance 
along the beam axis, or invariance to translations in spatial rapidity, is assumed. Spatial rapidity is a measure of the 
position along the beam, or z, axis, z/t = tanhr?. If all fluid elements begin at z = t = 0 and if the elements do not 
accelerate longitudinally, which is expected for a boost-invariant system, the spatial rapidity can be associated with 
the longitudinal velocity of the fluid element, v z = z/t. If a particle moves with the fluid, its spatial rapidity stays 
constant j7]. 



3 



FIG. 1. Charge fluctuations from lattice gauge theory [5j [6] (open symbols) are similar to those of a hadronic gas (filled 
symbols for T = 165 MeV). For fixed entropy there are increased numbers of up and down quarks in the hadronic phase, 
whereas the number of strange quarks is slightly smaller. The off-diagonal element disappears above T c when hadrons dissolve 
and quark-antiquark correlations disappear. At high temperatures the results approach those of a Stefan-Boltzmann gas of 
massless partons (S.B. limit). 


If chemical equilibrium is attained and if the correlations are local the correlation becomes 

9ab( A»j) ->■ Xabd(Ag), (5) 

where the S function can be relaxed to some function of short range that integrates to unity. For the realistic situation 
where the charge is created locally and diffuses over a finite distance, the correlation becomes 

gab(Ari) = XabS(Arj) + g' ab (Ag), (6) 

J dAg g' ab {Ag) = -\ a b- 

The last condition derives from charge conservation. 

Even though the net charge over the entire collision volume does not fluctuate and even though the short-range 
correlation only carries information about the current value of x, one can gain insight into the temporal history of 
Xab by analyzing the dependence of Xab on A 77 , or in a three dimensional analysis on Ax, Ay , Ay. For example, if 
Xab/s were to stay constant after the initial time, g' ab (Ar 7 ) would be a broad function with a width determined by 
how far charge pairs produced early would separate over the history of the collision. This width might be driven by 
a combination of two effects. First, if the balancing quarks are produced by the fragmentation of a longitudinal flux 
tube, the tunneling would lead to the charges being significantly separated at birth along the beam direction, perhaps 
by a solid fraction of one unit of spatial rapidity. For example, if the quarks were born with a separation of 0.5 fm, and 
were created at a proper time of 0.5 fm/c, the separation would be on unit of spatial rapidity. In contrast, if particles 
were born 0.5 fm apart at a time of 5 fm/c, the separation would be only 0.1 units of spatial rapidity. If the system 
maintains local chemical equilibrium, as defined by the lattice calculations displayed in Fig. [TJ at hadronization a 
second contribution to g a b would arise to account for the change in Xab • This second contribution to g' ab ( A 77 ) would 
be more tightly constrained in Ay due to the reduced time available for the charges to diffuse away from one another. 

Diffusion represents another mechanism of separation. For a strongly interacting QGP the diffusion constant and 
the separation would tend to be smaller. From [ 8 ] one can estimate the diffusive separation with a simple analytic 
formula based on assuming cross sections stayed constant with temperature. The diffusive width to the balance 
function was 


d) = 4^1n(r/r 0 ), (3 = v t /(nT<j). (7) 

If the cross section of a = 10 mb = 1.0 fm 2 were used, and if the density scaled by the time r were nr = 5 fm -2 , 
and if t/to = 10, the particles would separate by approximately one unit of spatial rapidity. If several collisions were 
required to re-thermalize a particle the diffusive separation would increase. 

The correlations in terms of charges, g a &, translate into correlations between specific hadronic species a and /?„ 

G a f}{Arj) = ((n a ( 0 ) - n a ( 0 ))(n^(Ajy) - n ) §( Arj))). ( 8 ) 

A relationship between g^Ay) and G a p(Ay) can be derived if one assumes that a small balancing charge is spread 
randomly, or thermally, amongst the hadrons |§l. 
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Using a blastwave prescription, which provides a parametric description of thermal emission overlaid with collective 
flow and includes the effects of decays, one can then project the correlations in coordinate space to correlations in 
momentum space. These correlations are divided by the yield of particles of type a and are called charge balance 
functions, 


B a p(Ay) = 


((n„( 0 ) - n & (y))(np( 0 ) - n 0 (Ar]))) 


Via T ^la 


(9) 


which have been measured by experiment 

In the next section, the methods are reviewed for translating charge correlations, g a b(Ar]), to G a p(Arj) and then 
to balance functions B a p(Ay) using a blast-wave prescription. Preliminary measurements from STAR are reviewed 
in Sec. urn while Sec .s |IV| and |V| illustrate how individual balance functions are determined by specific parameters 
in the model. In Sec. |V| eight model parameters are systematically varied in a Markov Chain Monte Carlo (MCMC) 
exploration of parameter space that implements model emulator techniques. The results make a strong case that the 
matter created in central Au+Au collisions at RHIC is close to being chemically equilibrated early in the collision. 


In Sec. VI these findings are summarized and some ideas for future analysis are presented. 


II. GENERALIZED BALANCE FUNCTIONS AND THE TWO-SURGE BLAST-WAVE MODEL 


In this paper we present comparisons to a simple model that can reproduce the behavior expected from the trends 
seen in the lattice results displayed in Fig. [l] The picture assumes two surges of charge production: one where the 
initial equilibrated matter would be formed and a second corresponding to the jump in the susceptibilities near T c . 
To that end, the model will assume that the first surge results in charge correlations significantly spread in rapidity, 
while the second surge, occurring as the temperature falls below T ~ 200 MeV, would result in significantly more 
compact correlations. The spread in spatial rapidity at breakup from the two surges will be described by two widths, 
a a for the first surge and gb for the second surge. First, we repeat the derivations in [9 and describe how these 
correlations carry over into correlations in rapidity between two specific hadronic species. 

The 3x3 charge correlation matrix in coordinate space, g' ab (Ar /), determines the observable hadronic correlations 
in the final state if one assumes that the differential additional charges in a small volume are distributed thermally 
amongst the hadrons. Here, we define a correlation between two hadronic species a and (3 [9], 

G a0 (Ar]) = ((n Q (0) - n a (O))(n 0 (Ar]) - n 0 (Ar]))), (10) 

gab{Arj) = ^ V G «/3( A v)(laaq/3b, (H) 

a/3 

where q aa is the charge of type a carried by the hadron species a. There are many more hadronic species than charges, 
so G a p cannot be uniquely determined from g without making assumptions. The factor of 1/4 accounts for the 
double counting of species. For instance, the sum over a,/3 includes both 7 r + 7 r - and 7 r - 7 r + . 

If one assumes that an additional differential charge density Sp a within a given volume element is spread amongst 
the various species randomly, i.e. thermally, the change in each species’ yield can be found by finding the chemical 
potential Sp a required to generate 5p a . 

Sn a = (e 5 ^ aqoia — 1 )n a (12) 

3 PaQaa^a^ 


where Sp has absorbed the 1/T factor typically used in defining chemical potentials and is dimensionless, 
refers to the average particle density. The three values Sp a can now be found from the three constraints in 
In turn these then determine Sn a . 


while n a 
Eq. (flO). 


3Pa — ^ ^ Qaa ^ ^ ^aQab^ Pb-> 
a b 

= VUb ad) ^6, 

b 

b 

5n a = fig Y. q a aX ( ab’ d) ~ ls Pb■ 

ab 


( 13 ) 
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Using the fact that in Eq. (10) n a — = 2 5n a , 

G a fj(Ari) = 4 «a?aaxi6 ad)_1 56c(A7?)Xcd ad)_1 a/3<i%- 


abed 


Using the expression for ^( had ) in Eq. (jsjl one can see that the form is consistent, 

j ^ ' Gg/3(Ari)q aa qi3i) = g a b{Ari'). 


(14) 


(15) 


a/3 


In the two-surge model, calculating the experimentally measurable G a p(Ay), where Ay is the correlation in relative 
asymptotic rapidity, involves the following steps: 


1. Calculate g^Ag) at hadronization. Assume a form consisting of two surges, 

e -Ari 2 /2a 2 B 




(16) 


(QGP) e Arl 
~ Xab (:27072)1/2 


Here, cr a describes the width in spatial rapidity to which charges created in the first wave have separated by the 
time one reaches breakup and gb characterizes the separation from the second wave at breakup. The strength 
of the first surge, X^ GP ^, is the susceptibility using a number per unit rapidity, or taken from lattice, 


(QGP) = X £ ttice) (r » 300 MeV) dS 
Xab s dg ’ 


(17) 


and xi h 6 ad) is given by the properties at chemical freezeout, which for this calculation will be defined by a tern- 
perature T c h e m = 165 MeV. One can also calculate the entropy density of a hadron gas at chemical equilibrium, 

^chem • 


2. G a/ 3 (Ar /) is calculated from g^ using Eq. (14). After substituting Eq. (16) into Eq. (14), 

e -Ar] 2 /2a 2 A e ~Ar] 2 /2a 2 B 

Gap^AT)^ — WA : a/3^a^l3 YTy^. 2 M/2 ^ B ,a(3B a Tl/3 


{2na 2 A ) 1/2 

W A ,a0 = ( X (had)-l x (QGP) x (had)-l^ 


ab 


ab 


(27rcr|) 1 /2 ’ 

Q/Sbj 


(18) 


W B 


, \ ^ (had) —1 

,a 0 = 4 2_ daaXab “ WA,a/3- 


ab 


3. Given G a p(Ag) at freezeout, the correlations in coordinate space are then mapped to correlations in spatial 
rapidity, G a p(Ay ), according to the blast-wave prescription. This mapping is done by the following Monte 
Carlo procedure, (a) Pairs of particles are chosen with probability proportional to n a rip. (b) The particles 
are then placed in coordinate space with the relative position randomly according to a Gaussian width a a • (c) 
The particle’s local momenta are generated thermally, then boosted transversely according to the blast-wave 
prescription described below, (d) If the particles are both stable the balance functions for the given relative 
rapidity bin is incremented by WA,apntot e a e f3> where n tot is the number of hadrons per unity rapidity and e a is 
the efficiency with which particle a is measured, and depends on the particle type, its transverse momentum 
and rapidity |2T (e) If one or both particles are unstable, they are decayed and the relative rapidity is calculated 
for each of the products a' and /?'. Each correlation G a ’p' is then incremented by WA,apntot e a' e (3' • Steps b-f are 
repeated using the weight wb and the width gb • Decays are accounted for by considering the intra-correlations 
between any of the particles of type A and species a' and a" and using a weight n a . 

4. After using the previous steps to calculate G a p(Ay), one can then calculate the “generalized” balance function, 
i.e. the balance function for specific species, by the relation 


B&{3{Ay) = 


Gqg(Ay) 

V U/Q; 


(19) 
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The denominator requires generating particles with the species being picked proportional to n a . Then one 
should boost and decay the particles, and increment the denominator for the type a' by the efficiency for the 
final species by an amount e a /7i( tot ). As long as the number of samplings is the same as the previous step, the 
normalization should be correct. 

One can test calculations along the way by considering the normalization of various quantities with perfect 
acceptance. 


^2(wA,a0 + WB,ap)npq0a = 4 q aa , 

0 

X / dAr] q aa G a0 (Ar])q 0b = -4\ a b, 

a0 

X] / dAy B a p(5y)q 0 b = -2 q ab , 

0 J 


( 20 ) 


( 21 ) 


The quantities also have the symmetry, G a p = —G a p and B a p = —B a p. 

The blast wave used to transversely boost the particles uses four parameters. The first is the kinetic temperature at 
breakup. Although particle yields are determined by the temperature T c h e m = 165 MeV, where chemical equilibrium 
would be lost, the particle’s momenta are determined by the kinetic temperature T^ n = 102 MeV. The second 
parameter is the transverse collective velocity, u_ j_ = 0.732. After being generated thermally, they are boosted by a 
transverse velocity chosen from the distribution, 


dN 

du x du y 


o ( U x~^ U y)/^ lU ± 


( 22 ) 


The two parameters u± and T^ n are chosen to roughly reproduce the mean transverse momentum of pions and protons 
seen by the STAR [22j and PHENIX [23] collaborations at RHIC. The third and fourth parameters determine the 
chemistry. In addition to T c hem, a parameter F B = 2/3 is used to reduce the yield of all baryons to account for 
baryon annihilation below T c h e m EH 25,- Rather than employing the rather large baryon annihilation factors here, 
calculations were also performed using lower values of T c h em as suggested in HI EH |27]. If the final p/it ratios are 
similar, results change little. 

The blast wave picture used here is perhaps the simplest model that can readily incorporate two surges while 
providing a reasonable mapping between the correlations in coordinate space and those in asymptotic relative rapidity. 
If both widths a a and a B are set equal to one another, it becomes a one-surge model. Of course, both charge creation 
and emission are more complicated than what can be described in this picture. First, although one expects two 
surges, the susceptibility seen in Fig. [l] does rise slowly as the temperature falls in the region T > 200 MeV. One 
could replace the two surges with a continuum of creation proportional to d/dr(x a b/ s )- Each point in space would 
then contribute to the final g a b- However it is clear that, if the system stays close to chemical equilibrium, the bulk 
of these contributions comes from two surges. Additionally, if either surge is not strongly confined to a specific time, 
the shape of the ensuing correlation should no longer be Gaussian, even if the separation is purely diffusive. A second 
weakness of the model comes from the assumption that the matter disassociates simultaneously according to a single 
breakup temperature, Tki n , but with yields determined by a single temperature T c h e m- For shorter lived resonances like 
the p meson, yields should substantially drop between the two temperatures. Thus, the contributions from resonance 
decays are probably significantly overstated. However, if a neutral particle decays and rethermalizes, the ensuing 
relative rapidity of the two charges is not wholly different than if the two decay prodcucts escaped unscathed, because 
these decays have relative momenta of typical thermal scales. Despite these shortcomings, the two-surge model should 
provide an insightful point of comparison. A more realistic picture of creation, transport and emission of conserved 
charges is currently being pursued by some of the authors. 


III. SUMMARIZING EXPERIMENTAL MEASUREMENTS 

The STAR Collaboration has presented preliminary balance functions for four combinations of species: 7r + 7r - , 
K + K~, pp and K~p [28, 29j, displayed in Fig. [ 2 ] Before discussing the result, we apply an approximate acceptance 
and efficiency correction so that the true widths of the balance functions can be better discerned. This correction 
is approximate and based on an assumption that the balance functions depend only on the relative rapidity, or 
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equivalently that the balance functions do not depend on the transverse momenta of the two particles. Assuming 
equal numbers of particles and antiparticles, the balance function can be exactly expressed as 


B(Ay) 


[dPadygdfb^ 6(1 Vb - Va\ ~ Ay) [N a p(p a , yaPPWyp ~ N a 0 (pa,y a ,pb,yb)] 

J dp a dy a N a (p a ,y a ) 


(23) 


with p a referring to the transverse momentum coordinates of the particle of type a and pb doing the same for type {3. 
The rapidities of the two types are then y a and y^. In terms of the balance function, 


N a p(p a ,y a ,Pb,yb) ~ N a p(p a ,y a ,Pb,yb ) = N a (p a ,y a )B (v \p a ,y a ,p b ,y b )Ay(p b ,y b ), 


(24) 


where Ap(pb,yb) is the efficiency that can vary between zero and unity, and is the balance function if the 

acceptance for the particle of type (3 were perfect. Now, to make the approximation, one assumes that an addi¬ 
tional particle of type /? at y 5 has its transverse momentum distribution distributed according to the single-particle 
probability, Np(pb, yb), with half the strength coming at yb = y a + A y and the other at yb = y a — Ai/. 


B^\p a ,y a ,p b ,y b ) « \b^{A V ) 

2 J dp b Np (p b ,y b ) 


(25) 


One can now insert Eq.s (24) and (25) into Eq. (23) and see that B^ p \Ay) factors out of the expression, 


B(Av) - r( p)(Av) I / dPadygNaipa, y_MUb = Va + A y) 
[ V) 1 yj \ 2 f dp a dy a N a (p a , y a ) 


(26) 


| [dpa/paNjj^ a,ya)A(y b = Va ~ A ?/) 1 


2 J' dp a dy a N a (Pa ' Va) 


Ap(y b ) = 


fdp b N { P(p b ,y b )A 0 (p b ,y b ) 


f dp b Njp(p b ,y b ) 


If the efficiency Ap(pb,yb) is perfect, the balance function is B^ p \Ay). The expression in the brackets on the r.h.s. 
represents the acceptance correction, C(Ay). 


B^(Ay) 

C(Ay) 


BjAy) 

C(Ay) ’ 

f dpadygNaiPa , ygjAyjyg + Ay) 
2 / dp a dy a N a (p a ,y a ) 


J dpadyaN a (pa,y a )Ap(y a - Ay) 
2 / dp a dy a N a (p a ,y a ) 


(27) 


The acceptance probability Ap(yb) is basically the ratio of the measured yield of particles of type j3 at yb relative to 
the true yield. Using a STAR acceptance and efficiency filter [21] . this can be found by generating particles with the 
blast-wave prescription, then seeing what fraction are recorded. The correction factor C(Ay) then requires averaging 
Ap over the various possibilities for y a , which can again be performed with the blast-wave prescription. Figure [ 2 ] 
shows both the original preliminary balance function from STAR and the corrected version. For large Ay the factor 
is large and the experimental uncertainties are magnified. Points are not plotted where the uncertainty surpasses the 
maximum size of the balance function. 













0.80 





0.0 0.5 1.0 1.5 


Ay 

FIG. 2. Preliminary balance functions from STAR (green 
circles) are shown from four pairs of species: 7r + 7r _ , K + K ~, 
pp and K~p. The upper curves (red squares) show the re¬ 
sult after applying the approximate acceptance and efficiency 
correction described here. Measurements are for the most 
central collisions, 0-5%. 




Ay 


FIG. 3. The acceptance and efficiency corrected balance 
funtion, B a p(Ay), from model calculations (red circles) is 
compared to the corresponding calculations with the accep¬ 
tance and efficiency for particle type (3 being perfect (blue 
triangles) and for both particles having perfect acceptance 
(green squares). If the assumptions, upon which the accep¬ 
tance and efficiency corrections were based, are justified, the 
three calculations would be nearly identical. 


Because the form for the acceptance and efficiency corrections was built on what may be a dubious assumption, 
they were also applied to model balance functions. In this case, the corrected balance functions can be compared to 
those calculated with perfect acceptance to see whether the difference between corrected and perfect balance functions 
differ. Figure [3] shows three sets of balance functions B a p(Ay), the corrected version as described above, the balance 
function where the efficiency for the particle of type (3 is perfect but where the acceptance and efficiency for a is 
given by the filter, and finally a balance function where the acceptance for both types is perfect. Because a balance 
function’s denominator divides out the efficiency of the first particle, all three would be identical if the balance function 
were truly only a function of Ay as assumed in the approximation described by Eq. (25). The parameters used in 
the model were for the two-surge model, chosen to roughly reproduce the experimental balance functions: a a = 1-0, 
gb — 1/3. The susceptibility for the quark gluon plasma was chosen to be ^2 a Xab/s = 0.18 and Xss/Xuu = 0.93, to 
be roughly consistent with the lattice calculations in Fig. [l]for temperatures near 300 MeV. 

Even without a model the STAR measurements make it clear that describing the separation of balancing charges 
with a single scale in relative spatial rapidity cannot describe the data. The four experimental balance functions 
for the most central collisions are displayed in Fig. |2j In the absence of acceptance corrections, a one-surge model, 
such as Eq. (18) with a a = ctb , would give balance functions of equal width if they were calculated as a function 
of spatial rapidity. However, given that the mapping to spatial rapidity is thermally smeared to a larger extent for 
lighter particles due to their higher thermal velocity, one would expect the following hierarchy of widths, a a p\ 


^TT+TT- > &K+K- > VpK- > Vpp- ( 28 ) 

The STAR data violates this expectation as it appears that 

Gpp > g k + k - > g^+t,- > G pK -. (29) 

The 7T + 7r - balance function should be strongly affected by decays. If one were to decay a chemically equilibrated 
hadron gas at T = 165 without any further hadronic interactions, nearly half the positive pions would come from 
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Ay 


Ay 


FIG. 4. Balance functions for the one-surge blast-wave model characterized by a width a a = 0.25 (dotted blue lines), a a = 0.5 
(dashed green line) and a a = 1.0 (solid red line) are compared to preliminary STAR results in the left-side panels. A width of 
~ 1.0 is necessary to explain the K + K~ and pp balance functions, but this width overestimates the width of the tv + ty~ balance 
functions and fails to approach the K~p balance function. The calculations include decays that occur below T = 165 and result 
in a narrow contribution that is similar but smaller than what appears in the measured 7r + 7r - balance function. The right-side 
panels show results for a two-surge model with a a = 1.0 and o B — 0.4. The two contributions to the K~p balance function 
have opposite signs which allows the resulting sum to be narrower than either the pp balance functions. The strength of the 
first surge is determined by matching the expected quark susceptibility to entropy ratio as extracted from lattice contributions. 


decays where the 7r + is accompanied by a 7r - . Thus, it should not be surprising that the 7r + 7r - balance function has 
a relatively narrow component. This could either be considered a third surge, or as an extension of the second surge. 
Determining what part of this narrowing is due to decays vs. the rise in the charge susceptibility in the transition 
region requires detailed analysis. However, the fact that the K~p balance function is narrower than either the pp or 
K + K~ balance functions, and even narrower than the 7r + 7r - balance functions, seems impossible to describe with a 
one-surge model. 


IV. COMPARING THE TWO-SURGE MODEL TO PRELIMINARY STAR RESULTS 


The left-side panels of Fig. [4] compare calculations using one scale, cta = cr B , at three values, a a = 0.25, 0.5 and 
1.0. For a a < 1.0 the model predictions provide poor descriptions of the data. For a a = 1.0, the K + K~ balance 
function is well produced and the pp balance function only differs at smaller relative rapidity. This latter failure could 
well be due to the lack of baryon-annihilation in the one-surge picture. The 7r + 7r - balance function from the model is 
wider than the data, as it seems the model would fit better with the narrowest of the three sample balance functions. 
As expected, the K~p balance function cannot be explained with a single scale. As expected, the single-scale model 
predicts widths that fall between the widths for the pp and K + K~ balance functions. Resonances that decay into 
both a proton and kaon are already considered in this calculation. 
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0.0 0.5 1.0 


Ay 

FIG. 5. The sensitivity of balance functions to various pa¬ 
rameters in the two-surge blast wave model is shown above. 
The default calculation (solid red line) assumes the first surge 
achieved chemical equilibrium according to lattice calcula¬ 
tions for T rsj 300 MeV. This would be a strangeness to up 
quark ratio, Xss/Xuu = 0.93, and a net quark to entropy 
ratio, (Xuu + Xdd + Xss)/s = 0.18. In the default calcula¬ 
tion the spreads of the two surges are cr a = 1.0, crs = 0.4 
units of spatial rapidity. The sensitivity of the K + K~ bal¬ 
ance function to changes in the strangeness content of the 
first wave is illustrated by comparing to calculations with 
Xss/Xuu = 0.5 (green dashed line) and Xss/Xuu = 1.4 (dot¬ 
ted blue line). For pp, calculations were performed for dif¬ 
ferent quark contents, (xuu A Xdd A Xss) — 0.1 (green dashed 
line) and (xuu + Xdd A Xss)/s = 0.25 (blue dotted line). 
K~p balance functions results are shown for varied widths, 
gb — 0.1 (green dashed line) and gb — 0.7 (blue dotted 
line). 
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FIG. 6. Balance functions for the two-surge model with de¬ 
fault parameters (solid black line) are divided into three com¬ 
ponents. The first surge (dashed red line) gives the widest 
contribution, with the width determined by a a plus some 
thermal smearing. The width from the contribution from 
the second surge (dotted green line) is set by as, plus ther¬ 
mal smearing. The third contribution is due to decays (solid 
blue line). The decay contribution provides nearly half the 
7r + 7r - balance function, but is much less important for the 
other species pairs. The K + K~ balance function is domi¬ 
nated by the first surge. The pp balance function has contri¬ 
butions from both, with the small negative contribution from 
the second surge resulting in a very flat balance function at 
small Ay. The K~p balance function has contributions of 
similar strength but opposite sign. This results in a narrow 
peak at small Ay that is much narrower than either the pp 
or K + K~ balance functions. Further, at large Ay the K~p 
balance turns slightly negative. 


The two-surge model provides a far superior description of the data than a one-surge model. To illustrate the 
sensitivity of the predictions to specific parameters, Fig. [5] shows how some of the balance functions react to changes 
of the default parameterization (a a = 1,<Tb = 0.4, (xuu + Xdd A x ss )/s = 0.18, Xss/Xuu = 0.93). In the default 
parameterization the quark to entropy ratios from the first surge are set to be consistent with lattice calculations. 
Changing the ratio of strange to up quarks, Xss/Xuu, from the first surge significantly affects the K + K~ balance 
function as seen in Fig. [5] For smaller strangeness in the first surge, the second surge must then pick up the 
difference, which gives a larger peak at small Ay. 

Adjusting the net-quark to entropy ratio, ( Xuu + Xdd A Xss)/ s, of the first surge affects all the balance functions, 
especially the pp balance function as shown in the middle panel of Fig. [5] If the quark content of the first surge 
is below equilibrium, the second surge would have to compensate for it and the pp balance function would have a 
positive narrow contribution, which is not seen in the data. 

Of the four balance functions studied here, the K~p balance function is the most sensitive to the separation of scales, 
c tb/&a , and most strongly illustrates the need for two surges. This derives from the cancellation of the opposite-sign 
weights for the A and B contributions in Eq. (18). If a a were to equal gb the balance functions would nearly cancel 
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parameter 

cta 

<?B 

(Xuu ± Xdd ± Xss)/ S 

Xss /Xuu 

Tkin 

u± 

Fb 

Avisc 

Min 

0.3 

0 

0.05 

0.0 

75 

0.5 

0.6 

0.7 

Max 

1.5 

1.0 

0.35 

1.3 

120 

0.875 

0.8 

1.0 


TABLE I. Eight parameters were simultaneously varied between the Min and Max ranges shown above. The first four 
parameters were two spatial rapidity widths, a a and (Jb, the quark to entropy ratio of the first surge (xuu + Xdd + Xss)/s , 
and the strangeness to up ratio of the first surge Xss/Xuu • The other four parameters vari ed the blast wave kinematics: the 
kinetic freezeout temperature Tki n , the transverse flow parameter u± defined in Eq. (22), the baryon reduction factor Fb 
which accounts for baryon annihilation, and A v isc, which accounts for the anisotropic shape of the phase space distribution in 
momentum space due to viscous effects. 


as can be seen in Fig. [4] As (Jb/va 0, both the narrow positive peak near Ay = 0, and the broad negative feature 
for Ay ~ 1 become more pronounced. 

The 7r + 7r - balance function is the least sensitive to changing model parameters for two reasons. First, the pions 
are light and their thermal motion smears any structure in relative rapidity. Second, nearly half the balancing 7r + 7r - 
pairs come from resonance decays that produce both a positive and negative pion. Figure [6] shows the contribution to 
each balance function in the default two-surge calculation from the first and second surges, and from decays below the 
chemical freezeout temperature. The 7r + 7r - balance function is strongly affected by decays. The decay contribution 
is constant in the context of this model, but could easily be sensitive to alterations of the model not considered here. 
In particular, this model assumes the chemical yields of various resonances do not change below chemical freezeout, 
T = 165 MeV, and kinetic freezeout when the particles are emitted, T = 102 MeV. In a realistic calculation resonances 
like the p would decay and their products would rescatter before final emission. Their contribution might then move 
to lower values of Ay because the relative momentum of the pions from a po decay tend to be higher than those 
from thermal motion. Of course, the yields themselves might be different. An improved calculation would include the 
effects of rescattering and baryon annihilation, better account for the spectral shape of the broader mesons, and more 
accurately describe the acceptance and efficiency for which weak decays, e.g. the K s , are captured by the detector. 

Because the strangeness susceptibility is rather flat near the transition region, and because few resonances decay 
to K + K ~, the K + K~ balance function is dominated by the first surge. In both the preliminary STAR data and 
in the models, contributions from 0 decays have been removed. The K + K~ balance function is then an excellent 
candidate to determine cr^. Due to baryon annihilation near and below chemical freezeout, the baryon susceptibility 
drops and the second-surge contribution to the pp balance function is negative. If not for thermal smearing, one would 
see a dip in the pp balance function. In reality, much of the annihilation occurs at the very end of the reaction and 
involves baryons with small relative momenta. Thus, the negative contribution from baryon annihilation should be 
more pronounced than the calculations shown here. The K~p balance function is especially insightful because the two 
contributions have similar strengths and opposite signs. For this reason, the resulting peak can be remarkably narrow, 
and followed by a shallow negative contribution at larger Ay. This makes it possible for modest changes in model 
parameters to significantly affect the size and shape of the K~p balance function. The negative dip at larger Ay 
requires a large experimental acceptance. Even though STAR has an acceptance of ±0.9 units of pseudo-rapidity, the 
acceptance in rapidity for protons and kaons is effectively narrower due to the mapping of rapidity to pseudo-rapidity 
for more massive particles. Combined with the small size of the balance functions, STAR’s measurements of the K~p 
balance function for Ay > 1 is difficult. 


V. STATISTICAL COMPARISON OF MODEL TO DATA 

Figure [5] makes a case that the first surge of particle production is close to what one would expect from chemically 
equilibrated matter. However, the study of model responses shown in Fig. [5] is quite incomplete given that model 
parameters can be changed in a coordinated fashion. A more rigorous method is to simultaneously vary all model 
parameters while comparing to all the data. Here, we present results from not only changing the four model parameters 
mentioned before, but four additional parameters. The posterior likelihood is calculated in this 8-dimensional space 
by comparing to all four balance functions. This will provide a distribution of likely parameters that not only provides 
the most likely point, but provides the full likelihood in the 8-dimensional space by a Markov Chain Monte Carlo 
(MCMC) procedure. 

Table [I] lists the parameters varied in this analysis. The first seven parameters were described in Sec. [IT] and the 
eighth accounts for the fact that the momentum distribution might not be anisotropic. After generating thermal 
particles, the momenta are scaled by a matrix 


- 'V A XX A yy A/2. 


Pi=Pi + KjPj, Kz 


(30) 
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This factor reduces the thermal smearing when mapping the spatial and momentum rapidities and gives narrower 
balance functions. However, it does not affect the width of the decay contribution. 

Even though the calculation of the balance functions requires only ~ 10 minutes of CPU time, hundreds of thousands 
of points need to be calculated during the MCMC trace. Using the techniques of model emulators in m, the likelihood 
is calculated by interpolating principal components of the data from 1536 full model runs. The full model runs were 
performed at points spread throughout the parameter space according to latin hyper-cube sampling. These techniques 
are described in [sung. 

The statistical analysis assumes Gaussian forms for the likelihood, 

I>i model) - ya xp) )^ab(yl model) - irt/ 2 } , (31) 

a ) 

^ab 

where the experimental, ?/i exp \ and model, 7/i model ^(x), are compared relative to the uncertainty a a . The eight 
parameters are considered to have uniform prior distributions as listed in Table [I] The uncertainties a a encapsulate 
not only experimental uncertainties, but any uncertainties one might expect from missing physics in the model. This 
can be thought of as a systematic theory error. For a schematic model such used here, assigning uncertainties is ad 
hoc, which means the final likelihood distribution is also suspect. However, even with that caveat, the statistical 
analysis helps determine what parameters best fit the data and can also assist with understanding which observables 
are best constraining specific parameters of interest. For this analysis, each individual point in the balance function 
was used, except for the first bin in relative rapidity. This bin carries both experimental difficulties due to two-track 
resolution and is strongly affected by femtoscopic correlations, which are ignored by the model. The uncertainty for 
each point was defined as 


exp - 


Ca = \/ (0-12y a ) 2 + (0.001) 2 (32) 

effectively a twelve percent error. This seems rather large, but because neighboring points in a balance function provide 
similar information and are similarly susceptible to theoretical approximations, groups of points can be redundant. 
For example, if the same quantity is measured four times, a 12 percent error falls to a six percent error. So if each 
balance function were divided up into groups of four points, the error would translate into half this uncertainty. The 
MCMC method described here was repeated with several forms for the uncertainty. For smaller uncertainties the 
likelihoods became more compact, but the position of the maximum likelihood did not change significantly. In addition 
to the balance functions, the mean transverse momenta for pions, kaons and protons were treated as observables. The 
values were averaged from STAR and PHENIX measurements and a 6 percent uncertainty was assigned. Each balance 
function had 17 points, due to the binning being in units of 0.1 in relative rapidity and extending to Ay <1.8, with the 
first bin being neglected. The four balance functions thus provided 68 observables, adding the three mean transverse 
momenta result in 71 observables. 

Figure [7] displays projections of the likelihood distribution as taken from the MCMC trace. The eight plots along 
the diagonal show one-dimensional projections, while the off-diagonal plots show two-dimensional projections. The 
quark content of the first surge is close to expectations from the lattice calculations shown in Fig. [lj but might be 5% 
high. Given that entropy production during the collision might be at the ten percent level, the entropy used in the 
denominator here might be ten percent higher than the entropy at early times. Thus, it would be not surprising if 
the x/ s ratios extracted here would be 10 percent lower than the expected values for lattice calculations. This implies 
that the (xuu + Xdd + Xss)/s ratio may actually be ^ 15% higher than expectations and that the strangeness content, 
may be about 15% lower than such expectations. This variation is at the 1 — g level in the comparison here, but 
because defining the uncertainty was rather arbitrary, stating a discrepancy between an extracted value and a lattice 
value is correspondingly arbitrary. 

The widths of the two surges are strongly constrained by the analysis. The first surge’s width is close to unity, which 
is in line with expectations. The characteristic width of the second surge appears to be just below half the broader 
width, ctb/cta < 0.5. This differs from the previous analysis of [33], which found a preference for gb/cta < 0.3. This 
change is being driven by the inclusion of the K~p balance function. If neglected, there is a preference to fit the 7r + 7r - 
balance function by lowering gb • Thus, the two balance functions are somewhat at odds with one another, which 
suggests that some physics is being poorly described or missing. A possible culprit would be the way in which decays 
are treated. If the number of decaying particles is significantly reduced between chemical and kinetic freezeout, the 
7 T + 7r - balance function, which is strongly affected by decays, might narrow. The balance function with gb/cta ~ 0.45, 
would then be able to better reproduce both balance functions. 

As a test of the statistical procedure we consider an assortment of balance functions calculated from random points 
in parameter space. In the left-side panels of Fig. [8] the parameters are taken randomly from the prior, i.e. they 
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FIG. 7. One- and two-dimensional projections of the posterior likelihood for the eight-parameter two-surge blastwave model, 
as computed using preliminary balance functions from STAR. Red, green and blue lines represent 1 — cr, 2 — a and 3 — a 
contours. The gray lines show lattice expectations for Xss/Xuu and (Xuu + Xdd + Xss)/s. The quark content of the first surge, 
(Xuu + Xdd + Xss)/s appears about 10 percent higher than the value expected from lattice and the strangeness to up ratio is 
about 20% low. 


are random numbers generated between the minimum and maximum values listed in Table [T| On the calculations 
for the right side of Fig. [8] the parameters were taken randomly from the posterior distribution, i.e. weighted by 
the likelihood. These points were extracted from taking 40 points far away from one another in the MCMC trace. 
Whereas the balance functions on the left-side panels vary widely and often stray far from the data, the balance 
functions from experimentally constrained parameters on the right side come close to the data. The 7r + 7r - balance 
function varies only modestly as parameters are varied throughout the prior, which means it has limited resolving 
power, although this could change if the 7r + 7r - balance function turns out to be significantly sensitive to missing 
or poorly represented physics, e.g. reabsorption of hadronic decays into the medium. In contrast, the other three 
balance functions vary widely throughout the prior, as was expected from the studies presented in Fig. [5] These 
measurements then have much higher potential to discriminate between different sets of parameters. The reasonably 
good model-to-data match validate that the statistical method was effective in identifying the most likely regions of 
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FIG. 8. The left-side balance functions were calculated from 40 parameters randomly taken from the prior, i.e. they were 
randomly generated in the regions defined in Table [I] The parameters for the right-side balance functions were also randomly 
generated, but the choices were weighted with the experimental likelihood. This validates the statistical procedure and serves 
as a guide for how strictly the definition of uncertainties constrains the fits. 


parameter space. 

Given the large numbers of parameters and observables, it can be difficult to understand the degree to which given 
measurements contribute to constraining the parameters, or the degree to which a measurement, if changed by one 
sigma, would lead to new parameters. In [34] methods were derived to perform sensitivity analyses from the MCMC 
trace. One such measure is the quantity, d((xi))/dyji* p ^ keeping all other y^ a fixed. This addresses the question 
of how the average value of a specific parameter X{ in the posterior would change if a single observable y a were 
altered without repeating the entire MCMC with a new value of ?/i exp \ In the expression above the double brackets, 
((• • •)), refer to an average from the posterior, and single brackets (• • •) denote the average over the prior. Figure [9] 
shows the values of this partial derivative where the parameters Xi are scaled by the widths of their priors so that 
((x — (x)) 2 ) = 1, and the observables y a are scaled by their uncertainties cr a . The resulting derivative then describes 
how much the parameter would change as a fraction of its prior width when the parameter y a is increased by a a . 
Figure [9] presents a much more detailed set of information than what is shown in Fig. [5] 


VI. SUMMARY AND OUTLOOK 

Preliminary STAR measurements of charge balance functions for specific species provide strict constraints on any 
picture of the chemical evolution of the super-hadronic matter created in heavy ion collisions. In order to fit all four 
of the balance functions, the blast-wave model used here had to employ two separate surges of charge production. 
Balancing charges from the first surge appear to have separated by approximately one unit of spatial rapidity by the 
time the matter reaches chemical freezeout. Charges from the second surge appears to have spread a bit less than 
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FIG. 9. The ratio d((xi))/dy^ p ^ is shown for every combination of 71 observables and 8 parameters. This describes how 
an extracted value for a parameter, ((#)), would change if a given observable were increased by one sigma while keeping the 
other observables fixed. Each of the four balance functions provides 17 observables, corresponding to the 18 rapidity bins with 
the first bin being ignored. The first 17 observables along the x axis represent the 7r + 7r - balance function, and the other three 
balance functions are represented by each subsequent group of 17 observables. Additionally, the three mean transverse momenta 
for pions, kaons and protons are included as observables. For example, one can see that the extracted value of Xss/Xuu would 
fall/rise if the low momentum bins of the K + K~ balance function were raised/lowered. The response is opposite for bins of 
higher relative rapidity. One can also see the relatively strong resolving power for the K + K ~, pp and K~p balance functions. 


half that amount. It is difficult to determine the physics driving the spread of the first surge. It may mainly be due 
to the separation inherent to breaking strings, or equivalently the tunneling of balancing charges in the breaking of a 
flux tube. Another possibility is that the spread is largely diffusive. Either way, the charges must be created early if 
they are to separate by such a large distance by chemical freezeout. 

The strengths of the two surges can also be determined by comparing the two-surge model to data. Remarkably, 
the strength of the first surge is consistent with the matter from the first surge approaching chemical equilibrium. A 
model run based on assuming perfect equilibrium, followed by an isentropic expansion, reproduces the preliminary 
STAR data at the 5% level. A more comprehensive search through eight-dimensional parameter space, including four 
parameters that adjust the blast-wave description, show that the best chemistry for fitting the data has susceptibilities 
within 20% of those calculated on the lattice. The response of the measurements to the susceptibilities was found to 
be strong, and thus should provide good resolving power within the context of any model. Measurements were shown 
to be sensitive to both the spread and the strength of the charges. 

The preliminary STAR measurements considered here are the first of their kind, and should inspire numerous 
additional measurements of different species, measurements as functions of relative azimuthal angle [35, 36 and 
transverse momentum, and as a function of the three-dimensional pair momentum. As a general statement, balance 
functions are six-dimensional correlations that can be measured between any two species pairs. Every new combination 
of species and every more differential binning of the six-dimensional phase space should provide new insights and better 
shed light onto the story of the chemical evolution of the reaction. In addition to higher statistics, these goals may 
also require expanding the acceptance of current experiments. Here, we provide a quick list of some measurements 
that would help illuminate an assortment issues: 

• Measurements of unidentified particles have been made as a function of both the relative and total azimuthal 
angle. However, one could consider the relation between the width in azimuthal angle and relative rapidity. 
Balancing charges that are created early can be preferentially selected by considering the balance function at 
large relative rapidity, and should also have broader balance functions in relative azimuthal angle. This might 
help distinguish the physical cause of the separation of early charge pairs, diffusion vs. flux tube breaking. 

• The STAR Collaboration’s acceptance for identified particles extends to ±0.9 units of pseudo-rapidity, but 
when translated into real rapidity is significantly narrower for heavier particles like kaons and protons. The 
ALICE detector at the LHC covers ±0.8 units of pseudo-rapidity. ATLAS and CMS cover much wider regions, 
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but without particle identification. The long-range correlations should better illuminate the early production 
charges. For instance, if a measurable fraction of pairs comes from the initial state, those balancing charges 
might be particularly well separated. The greater acceptance would also make it easier to perform the analysis 
mentioned in the preceding bullet. 

• In order to better understand decays, balance functions can be binned as a function of the invariant mass or 
momentum. Contributions from specific decays, e.g. K s —» 7r7r, can be identified by their peak. Contribution 
from broader resonances such as po 7r + 7r - may be more difficult to see, but nonetheless the comparison 
should be able to identify whether the populations of unstable resonances are poorly represented by this, or any 
other, model. 

• Measurements can be performed as a function of beam energy and centrality m and for different choices of 
colliding nuclei [13]. The narrowing of the balance function for unidentified particles may be related to the 
existence of two surges [8], but that conclusion needs to be investigated in detail alongside the analyses listed 
above. 

If these measurements are to lead to truly rigorous conclusions, the models must improve beyond the schematic blast 
wave model considered here. In order to justify Gaussian forms for the two surges, each surge must be created during 
a short time. More realistically, charges are created continuously according to the rate at which Xab/s is changing. 
The evolution of the charges and the creation of new charge pairs after chemical freezeout is an area that was poorly 
modeled here. By the time of kinetic freezeout many of the resonances will already have decayed, with their products 
being reabsorbed by the medium. Baryon annihilation, which was implemented here by reducing the baryon yield 
at chemical freezeout, in fact occurs mostly later. The corresponding dip should then be more concentrated at low 
relative momentum, a feature that is indeed seen in the data. These improvements cannot be naturally added to a 
schematic blast wave model, but they can be accounted for in microscopic simulations. Such calculations are entirely 
tractable and are being performed currently by this group. Once the models can encapsulate all the important physics, 
statistical analyses like those presented here can lead to much more convincing and rigorous conclusions. 
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